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Abstract 

Many biological networks have to filter out useful information from a vast excess of spurious interactions. 
We use computational evolution to predict design features of networks processing ligand categorization. The 
important problem of early immune response is considered as a case-study. Rounds of evolution with differ- 
ent constraints uncover elaborations of the same network motif we name "adaptive sorting". Corresponding 
network substructures can be identified in current models of immune recognition. Our work draws a deep 
analogy between immune recognition and biochemical adaptation. 

PACS numbers: 87.16.Xa, 87.18.Mp, 87.18.Tt, 05.10.-a 
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Information processing in biology often relies on complex out-of-equilibrium physical pro- 
cesses ensuring efficiency [IJ. The paradigmatic example is kinetic proofreading (KPR), first 
proposed to explain low spurious base-pair interactions during DNA replication KPR orig- 

inated in a context with comparable concentrations of correct and spurious substrates. If the spu- 
rious substrate has similar characteristics and is orders of magnitude higher in concentration than 
the correct one, alternative strategies are needed. 

An important instance of this problem is immune recognition by T cells. T cells constantly scan 
antigen presenting cells (AFC) in their environment, via the binding of their T cell receptors (TCR) 
to the presented pMHC ligands. T cells perform a sorting process based on interaction with self 
(non-agonist) or foreign (agonist) ligands at the surface of APCs: if foreign ligands are detected, 
then the immune response is triggered. Following the "life-time" dogma [4J, one of the main 
determinants for distinguishing self from foreign is the unbinding time of the pMHC ligand to 
TCR. Ligands up to a critical binding time of Tc ^r:::^ 3 s do not elicit response while foreign ligands 
bound for a longer time (r/ > Tc) do. Self ligands dissociate rapidly (typically for < 0.1 s). 

The sorting process is extremely sensitive: response is triggered in the presence of minute con- 
centrations of foreign ligands (of the order of 1-10 ligands per cell JHO). Sorting is specific: 
although foreign (r/) and critical ligands (Tc) have similar binding times, an arbitrary concen- 
tration of critical ligands does not elicit response These requirements are summarized on 
Fig. [TJ McKeithan [8J proposed that T cells harness the amplifying properties of KPR to solve the 
recognition problem between few foreign ligands and vastly numerous self ligands. However, this 
model does not account for sharp thresholding required for sensitivity and specificity as noticed in 
Q. Other control structures must exist. 

We use computational evolution [9J to ask the related "inverse problem" question: how can a 
network categorize sharply two ligands with similar affinity irrespective of their concentrations? 
We exhibit networks performing ligand recognition with the help of a new network module that we 
name "adaptive sorting" which we study analytically. We use extensive evolutionary simulations 
to show how this solution is improved to solve the related recognition problem of parallel sorting 
of foreign ligands within a sea of self ligands. We expect the principles presented here to have 
broader relevance for biological recognition systems where specific signals must be extracted from 
a high number of weak spurious interactions. 



2 



(a) APC 




(b) 



T-Cell 



APC 



Response 



T-Cell 



3 



No Response 



(c) 



Specificity 



3 




J Sensitivity 



<0,1 3 10 Binding Time (s) 

Critical ^ Foreign T-Cell 
"(t, = 3s) w"(Tf>10s) >^ Receptor 



FIG. 1. (color online) Schematic illustration of the problem setup, (a) Few foreign ligand (r/ > 10 s) 
trigger response, (b) Arbitrary large concentrations of critical agonist (tc = 3 s) ligands do not trigger 
response, (c) Idealization of the number of pMHC ligand required to trigger response as a function of 
pMHC-TCR binding time. Shaded region corresponds to conditions for which the immune response is 
triggered. 

Methods. - The algorithm we use to generate biochemical networks is essentially the same as 
that used in fTO], with a biochemical grammar adapted to the specific problem of ligand recogni- 
tion by (immune) cells. Following the model described in [|7]], we limit our possible interaction 
grammar to phosphorylations or dephosphorylations with rates linear in enzyme concentrations. 
We assume ligands bind TCR outside the cell, resulting in the activation of the internal part of 
the receptor (denoted by Co, see Fig. [T] (a), (b)). The algorithm then proceeds to add/remove 
kinases/phosphatases to evolve cascades of phospho-reactions downstream of Cq. We make the 
classical hypothesis underlying KPR models [|8jl that when a ligand dissociates from a receptor, 
the receptor's internal part gets quickly dephosphorylated. This assumption is consistent with the 
"kinetic segregation" mechanism ifTDl (see details in US). We assume that a single species in the 
network plays the role of the output of the system and triggers immune response in a binary way 
via a thresholding mechanism. The nature of the output is under selective pressure and can be 
changed by the algorithm. 



The goal here is to discriminate between two kinds of ligands with identical on-rate (denoted 
by tv) but different binding times: r/ = 10 s for foreign ligands and Tc = 3 s for critical ligands 
(we checked that our results do not depend on the specific choice of r/ and Tc as long as both are of 
the same order of magnitude). For pure KPR [HI, the concentration of the output is linear in ligand 
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concentration. Thus, as shown in Fig. [2] (a), hgands with similar binding times are distinguished 
by a thresholding mechanism only over a limited range of concentration, even for a large number 
of proofreading steps [7J. In contrast, if the steady state output concentration is almost flat in 
ligand concentration due to some control mechanism, as shown in Fig.[2](b), then ligands can be 
categorized by thresholding nearly irrespective of their concentration. 




FIG. 2. (color online) (a) KPR scheme has discrimination abilities over a limited range of ligand concen- 
tration, (b) Output ligand for rj = 10 s and Tc = 3 s. (c) Histogram of outputs from (b) illustrating 
effective probability distribution, (d) Adaptive sorting network. Arrows with no specified enzyme represent 
unregulated reactions. The output is circled. We keep conventions throughout, (e) Output ligand and 
histogram of output for adaptive sorting = 10~^, R = 10^, 5 = 1, e = 1, a = 0.3, 6 = and Kt = 1). 
(f) Minimum ligand concentration triggering response for different binding times for adaptive sorting in (e). 
Threshold taken to be C(^c)- 

To select for networks producing almost flat ligand dependency, we start by sampling loga- 
rithmically the range of allowed ligand concentration. Then, steady state outputs are computed 
for every ligand concentration and binned for the two binding times considered (Fig.|2](c) shows 
the binned outputs corresponding to Fig. [2](b)). One then considers the histograms of output for 
different r's as an effective probability distribution function. A natural measure of performance 
("fitness") selecting for networks with behaviour similar to Fig. |2] (b) is then the mutual infor- 
mation, X(0;r = {tc, r/}) [|T3ll , between the output value and the dissociation time. A network 
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for which X{0\t) = Imax{= 1 bit) has its output distributions for r/ and Tc disjoint, which is 
biologically equivalent to a perfect discrimination. We run our evolutionary simulations with this 
fitness function for 30 values of ligand concentrations equally spaced on a logarithmic scale in the 
interval [1 10^]. More details on the evolutionary simulations are given in lfT2l . 

Simple adaptive sorting. - We run our simulations with deterministic integration of network 
equations. Figure [2] (d) presents a typical network topology we obtain, with the corresponding 
distribution of outputs on Fig. |2] (e). Distributions corresponding to the two binding times are 
clearly separated. In this network, Co is phosphorylated once into Ci by kinase K. K is itself 
phosphorylated by Co, which makes it inactive. Here, Ci is the output. Calling R, L and Kt the 
total concentration of receptors, ligands and kinase respectively equations for this network are 

Co = nR'^^'L''^' - {aK + r'^) Co + bCi, (1) 
Ci = aKCo-{r-^ + b)Ci, (2) 
K = -5CoK + e{KT-K). (3) 

^free = _ ^^1=0 Ci and L^^^^ = L — J2i=o concentrations of free receptors and 

ligands. Assuming receptors are in excess (i?^^^^ c::^ i?), the steady state concentration of output 
variable Ci can be easily computed and we get Ci = ^r^^ where ^{r) = f^^^, = e5~^. 

For large L, Co oc L. In particular, as Co^ C^.Ci ^ ^i^)- It is also clear that even for small 
L, Ci will be a pure function of r independent from L if small enough. To discriminate between 
two ligands with binding times ri and r2, one then simply needs to assume response is activated 
for a Ci threshold value 6 G [C(^i) , ^(^^2)] • Figure |2](f) illustrates the range in ligand concentration 
leading to response with such thresholding (taking 6 = ^{tc)) process for the present network. 
The network shows both extremely good sensitivity and specificity (compare with Fig.[T](c)). 

This situation is reminiscent of biochemical adaptation, where one variable returns to the same 
steady state value irrespective of ligand concentration. Indeed, the motif displayed on Fig. |2] (c) 
implements an "incoherent feedforward loop" logic as observed in adaptive systems lfT0l[T4l[T5]l : 
Co feeds negatively into kinase K, and both Co and K feed positively into output Ci. The overall 
influence of Co (and of L) is a balance between two opposite effects which cancel out. However, 
one significant difference from classical adaptation is that the steady state concentration of Ci is 
now a function of the extra parameter r, the ligand dissociation time. Discrimination of ligands 
based on the value of the output becomes possible irrespective of the ligand concentrations. 

This process can be generalized to other adaptive networks based on ligand-receptor interaction, 
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as long as one kinetic parameter is ligand specific. For instance, ligand-receptor networks evolved 
in ifTOl can be modified to have a steady state concentration depending on ligand nature. Call / 
the input and R the (stable) receptor. If we assume that the complex C resulting from association 
/ and R is washed away with a time constant r/ depending on the nature of the input, then the 
simple adaptive system R = p — IR and C = IR — C/ti stabilizes to a steady state concentration 
C = pr/, which depends only on r/ irrespective of input value. Schematically, / plays same 
role as Co (proportional to ligand concentration in Eq. [T]) while R plays the role of K (inversely 
proportional to Co from Eq. [3] and buffering it in Eq. [2]) . We believe that this combination of 
biochemical adaptation with a kinetic parameter dependency to perform decision can potentially 
be observed in a wide variety of biochemical networks. We subsequently call it adaptive sorting. 

Parallel adaptive sorting. - Adaptive sorting by itself is efficient to discriminate independently 
critical from foreign, but its performance is degraded when cells are exposed at the same time 
to foreign ligands (concentration Lf) and a huge excess of self ligands (concentration L^), as 
illustrated in Fig. [3] (a). This phenomenon is not specific to the immune system and is called 
antagonism [7]. The reason is that the two different kinds of ligands are coupled through the 
common kinase used in the feedforward motif (dashed arrows in Fig.[3](b)). Precisely, denoting 
the complexes arising from the binding of foreign and self ligands by Ci and Di respectively, the 
total output concentration is 

which still tends to i{rf) at large Lf. We can neglect Di in the output because ^(r) oc r and 
so ^{ts) <C i{Tf). To reach the adaptive regime, we now have the requirement that Co ^ ^o- 
For large Lg, Dq :» Di and we have that L^o + — ^ tvRTs{l + tvRTs)~^Ls. Similarly, 
Co ^ KRTf{l + KRTf)-^Lf. Thus Ci ^ for 
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With hiRvf :» 1, ^ 10^ and hiRrg ^0.1, self ligands thus annihilate the sensitivity of the 
simple adaptive sorting motif. 

To solve this problem, we rerun evolutionary simulations with the constraint that discrimina- 
tion between r/ and Tc should happen even in the massive presence of self ligands (r^ = 0.05 s), 
as sketched in Fig. |3] (c). A representative result of this computational evolution is presented in 
Fig. [3] (d) and (e) for output and network topology respectively. The networks found look very 
similar to adaptive sorting, except that the incoherent feedforward module is sometimes imple- 



mented via activation of a phosphatase, instead of de-activation of a kinase lfT6ll . A full cascade 
of KPR also evolves. Notably, in all working networks there is an important difference with the 
previous case: activation of the enzyme in the adaptive sorting module is rewired downstream the 
first step of the KPR cascade (dashed circles in Fig.[3](e)). 




FIG. 3. (color online) (a) Effect of self ligands on the adaptive sorting module from Fig. [2](e), taking 
[Cn] + [Dn] as an output. Full lines: Lg = 0. Dashed line: Lg = 10^. Note the catastrophic effect of 
self ligands on sensitivity (quantified through AL). We compare r/ with Ls > to Tc with = as a 
worst case scenario, (b) Coupling (dashed arrows) between two different types of ligands through kinase 
K for adaptive sorting, (c) Schematic illustration of new constraint of parallel sorting. Squares represent 
self ligands (r^ = 0.05 s). (d) Example of evolved output v^. ligand relathionship with = (full) 
and Lg = 10^ (dashed). Loss in sensitivity is now small, (e) Schematic of network corresponding to (d). 
Complexes C^'s are understood to decay to R^^^^ and L^^^^ (same convention in Fig.|4]). Parameters are given 

in ma. 

This can be fully understood analytically by considering an idealized network such as the one 
in Fig. |4] (a) which is compared to the actual network implicated in immune response [|7l [17]] in 
Fig.|4](b). Our idealization consists in an adaptive sorting module with upstream and downstream 
steps of KPR (A^ steps in total, adaptive module activated by complex m, m + 2 < N, Fig.|4](a)). 
In such networks, assuming no dephosphorylation down the cascade (b = 0), the output takes the 
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form iri2l 



Cm + D 



N 



c 



N 



(6) 



where Cm = jJ^Cq and = iT^o^ with 7^ = (/)Ti{l + (/)Ti)~^ . cp denotes the default (unregu- 
lated) phosphorylation rate in the cascade. ^\t) is a function of r, and like before ^\rs) <C ^\Tf) 
so that we can neglect the contribution of at in the output. Even in the presence of many self 
ligands L^, we clearly have an output independent of L/ for Cq 7/ ^T^^o ( m = is simple 
adaptive sorting). Since (/) ^ tJ^ for a sensitive network lfT2l . 7s is small, thus any m > 1 
makes jj^j^ even smaller. This is in essence a weak proofreading process upstream the cascade 
ensuring that Cm ^ Dm so that the adaptive sorting module is only triggered by foreign ligands. 
As for simple adaptive sorting, we have that Cq (x Lf and Dq oc Lg although the prefactors differ 
lfT2l . In the end, Cn is a pure function of Tf for 

so that the r.h.s is small compared to Eq. [5]for m > 0. Self influence is consequently almost 
abolished. 
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FIG. 4. (color online) (a) Final network with categorization properties in presence of large concentrations 
of spurious substrates. The parallel (long dash), adaptive (fine dash) and KPR (dotted) modules are identi- 
fied. Star indicates the specific phosphorylation in adaptive sorting, (b) Network for immune recognition 
with corresponding features, from ||7l[l7|. Adaptive sorting is achieved via the activation of non specific 
phosphatase P*, assumed to be SHP-1. (c) Sample trajectories of A for different ligand concentrations L. 
Warm color for r/ and cold colors for Tc. Threshold 9 is identified by an horizontal line. Black curves 
are the analytic expressions {A(t)) ± ^A{t) El- (d) Fraction of trajectories having reached threshold for 
TV = 4, m = 2, = 10-^, R = 10^, 5 = 1, e = 0.5, (j) = 0.3, a = 0.0003, 6 = 0, = 1000, Ls = 0. 
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It must be emphasized that the more complex solutions displayed in Fig. [3] (e) and Fig. |4] (a) 
require more than one kinase or phosphatase: generic enzymes are shared by most of the proof- 
reading steps, while a specific enzyme accounts for the adaptive sorting module (star in Fig. |4] 
(a)). This is of biological importance since it is not clear that biochemistry would allow fine-tuned 
specificity to a single step in the cascade. Interestingly, alternative solutions also evolve where 
kinases and phosphatases are not specific to a given proofreading step [il2il . For these networks, 
discrimination is still possible, but loss of biochemical specificity degrades the adaptive proper- 
ties. Instead, one observes a non-monotonic behaviour, flattened out over the range of input ligand 
considered, so that adaptation is only approximated. For a complete analytic and experimental 
study of such a case, see ifTTl . 

Dealing with low numbers of molecules - Immune cells perform efficient sorting of different 
ligand types for as little as ^ 10 foreign ligands. A low number of molecules is potentially 
problematic because adaptive sorting shows a trade-off between specificity and sensitivity. In the 
simpler scheme ( Fig. |2] (c)) , perfect adaptation for L ^ occurs if 0, but output in the 

adaptive regime is Ci = ^(r) oc ^ so that discrimination becomes impossible. Increasing 

actually softens the constraint: KPR steps downstream the adaptive module (Fig. |4]) add a 
geometric dependency in r to Cn , so that Cn can have a strong dependency in r (specificity) 
even for low (sensitivity) lfT2ll . 

Another potential problem comes from fluctuations at low ligands. In the immune context, 
fully phosphorylated tails of receptors (corresponding to Cn in our model) themselves slowly 
phosphorylate abundant (> 10^ molecules) downstream targets such as ZAP70 and ERK, which 
saturates and triggers response only after a couple of minutes [7J. Following ifTTl . we check that 
coarse-graining this downstream cascade into a slow variable solves the fluctuation problem. We 
pose a variable A obeying A = ACn — T~^A. ^4 is a proxy for the abundant targets and can 
therefore be realistically assumed to be deterministic as long as A is large, so that the only A 
stochasticity comes from Cn- We assume thresholding is then made on the deterministic value of 
A, leading to a binary irreversible decision ifTSl . We take T = 60 s, as the response of T-cells 
occurs on the order of minutes fl\. 

Simulations of this process using Gillespie algorithm are presented in Fig. |4](c) and (d), with 
samples of trajectories and fraction of activated cells as a function of time. Results are in very good 
agreement with a simple linear noise approximation on CAr(see details and assumptions in lfT2ll - 
in particular, fluctuations of A decay as T"^/^). Ligands at Tc essentially never cross the threshold 
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for the considered time window, while for hgands at r/, almost all cells eventually respond for 
Lf > 5 (this number goes to < 10 in presence of self ligands lfT2ll \ Finally, the model's half 
population response time (Fig. |4] (d)) is consistent with experiments (Tl [121 113 and decreases 
down to less than one minute as increases. So, although we cannot exclude that other noise- 
resistance mechanism are possible lfT9l , adaptive sorting coupled to a slow downstream cascade 
has discrimination capabilities compatible with experimental data. 

Many biological systems have to filter out specific useful information from a vast excess of 
spurious interactions. We have evolved in silico networks categorizing ligands with very close 
biochemical properties irrespective of their concentrations. We have discovered a new functional 
unit, the adaptive sorting module, which can be rewired to solve a parallel sorting problem. Our 
final model is summarized in Fig.|4](a), along with corresponding network features of the immune 
system Fig.|4](b) [fTTl . Strikingly, the network of the immune system shares many similarities with 
our final solution. In our framework, immune recognition corresponds to an optimal solution with 
non-specific enzymes. We expect adaptive sorting to manifest itself through non linear (or even 
non monotonic) dependency of response on input concentration. This is observed in a wide range 
of networks, for instance in endocrine signalling [20J, but remains mechanistically unexplained. 
Adaptive sorting could lie at the core of such signalling processes as well as others. 

We thank Eric Siggia, Massimo Vergassola, Guillaume Voisinne and Gregoire Altan-Bonnet 
for useful discussions. JBL is supported by NSERC, PF by NSERC and HFSP. 
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I. GRAMMAR OF ALLOWED CHEMICAL REACTIONS 



The starting point of all our evolutionary simulations is the minimal network shown below. 



where L and R are the total concentrations of ligand and receptor respectively, i?^^^^ = i? — Co 
is the free receptor concentrations and L^^^^ = L — Co is the free ligand concentration. The circle 
denotes that Co is the output of the network. Biologically, we consider Co to be the "activated" 
intracellular section of the TCR bound to a pMHC. 

The grammar of the network allows for the complex Co to be phosphorylated: 



We model the phosphorylations and dephosphoarylations in the simplest possible scheme. The 
equations for this network are then 



Here, Co gets phosphorylated to Ci by kinase K. Each phosphorylation is always associated 
with a dephosphorylation (to avoid irreversible reactions). In this network, the dephosphorylation 
is catalysed by phosphatase P. Observe also that Ci is taken to decay to free ligands and receptors 
with rate r"^. This reaction amounts to a kinetic proofreading step. The decay of Ci directly to 
receptors and ligands is equivalent to demanding fast dephosphorylation of the internal section of 
the TCR relative to the binding time r (the importance of this assumption is discussed in more 
details in Sec. |V]). Note that the output is left under selective pressure (i.e. the output tag is not 
carried to Ci automatically). The grammar allows for this reaction (phosphorylation of C^) to 
occur an unlimited number of times. 

We assume that complexes in the cascade (the C's) are kinases. This means that C^ can phos- 
phorylate any kinase or phosphatase, except other complexes C^. To have the reaction grammar 
as unbiased as possible, we take the catalytic activity of any kinase (or phosphatase) to initially 




With corresponding equation 



Co = /^(i?-Co) (L-Co)-r-iC( 



K 




Co = {Rt -C,-C,) {Lt - Co - Ci) - {r'' + aK) Co + /5PCi , 
Ci = aKCo- (/3P + r-^)Ci. 
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be shared with its phosphorylated counterpart (e.g., if initially K phosphorylates Co and K gets 
phosphorylated to X*, then X* also phosphorylates Co). The algorithm allows for the removal 
of interactions in subsequent generations if that proves advantageous in terms of the chosen fit- 
ness. To take a concrete example, if in the previous reaction Co phosphorylates K with rate 5, the 
network becomes as below. 



k = -5C^K + eM{K - Kt). 

Observe how K + = Kt, and not just K, enters the equation of Co and Ci. This inheritance 
of catalytic activity also holds for complexes in the kinetic proofreading cascade (the C's). 

As mentionned before, reactions can also be removed by the algorithm. For example, the 
reaction of X* with Co could be removed. The output tag can also shift between different species. 
For instance, the output could become Ci. One then ends up with network 



This is just the adaptive sorting module. 

It must be made explicit that our grammar assumes that the catalysis have fast kinetics. 
This means that an enzyme concentration is unaffected by its catalytic activities. Moreover, 




With corresponding equations 



Co = ^ [R^ - Co - Ci) {Lt - Co - Ci) - [r'^ + a{K + K*)} Co + /5PCi 
Ci = a{K + K*)Co - + r-^) Ci, 




In which case the equations read 



Co = {Rt - Co - Ci) (Lt - Co - Ci) - (r"^ + aK) Co + /5PC 
Ci = aKC^- (/3P + r-^)Ci, 
k = -5CoK + eM{K - Kt). 
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the strongest non-linearities are quadratic. This assumption is mainly motivated by simplicity. 
Even with these simple ingredients, solutions with interesting features are found. 

Finally, the grammar allows for more than one kinase to phosphorylate a given reaction and 
similarly for phosphatases. Also, there is a possibility to add kinases and phosphatase (we typi- 
cally start evolution with two kinases and phosphatases present). In addition to having dynamic 
topologies in the algorithm, all kinetic parameters as well as concentrations of enzymes not on the 
C cascade can be modified within predetermined range. 

Initial conditions for the integration of the networks' equations are taken to be as biologically 
realistic as possible. Before cell-cell contact, there clearly no bound TCR and pMHC. We therefore 
take = Vi at t = (the antigen presenting cells contacts the T cell at time 0). 

When performing evolution with cells exposed simultaneously to self and foreign ligands, we 
make coupling between self and foreign complexes as shown in the case of simple adaptive sorting 
in Fig. 3 (b) of the paper. If the output is a member of the cascade of complexes, then it is taken to 
be the sum of complexes arising form self as well as foreign ligands. We run our simulations for 
30 values of ligand concentrations equally spaced on a logarithmic scale in the interval [1 10^]. 

II. DERIVATION OF ASYMPTOTIC OUTPUT CONCENTRATIONS 
A. Simple Adaptive Sorting 

We consider the following system of differential equations (shown schematically in Fig. 2 (d) 
of paper) 





D,^aKDo-{T;' + b)D,, 



{Co + Do)K + e {Kt - K) . 



(8) 



The Cs and the D's are respectively the agonist (foreign or critical) and the non-agonist (self) 
complexes. The output is Ci+Di. The following tables summarize the meaning of each variable. 
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Symbol 


Chemical Species (Concentration) 


r 


Agonist complex phosphorylated n times 


Dn 

lb 


Non-agonist complex phosphorylated n times 


Ls 

b 


Self ligands, binding time = 0.05 s (total) 




Critical ae^onist lie^ands bindine^ time = 3 s ftotal) 


Lf 

J 


Foreign ligands, binding time r/ = 10 s (total) 


R 


Receptor (total) 


K 


Kinase 


Kt 


Kinase (total) 



Symbol 


Kinetic Parameter 


Order 


a 


Complex phosphorylation rate 


2 


b 


Complex dephosphorylation rate 


1 




Self complex binding time 






Critical agonist complex binding time 






Foreign complex binding time 




5 


Kinase phosphorylation rate 


2 


e 


Kinase dephosphorylation rate 


1 



We suppose that the receptors are largely in excess, i.e., (i? — XlLo + A]) ^ R We are 



interested in the steady-state concentration (see Sec. |VII| ). We first consider the case where the 
are no self ligands, Lg = 0. The case of two ligand types is treated in sec. |IV A[ Our goal is to 
determine the behaviour of Ci (steady- state) as function of Lf. Unless otherwise specified, we 
take r/ = 10 s, Tc = 3 s and Tg = 0.05 s. 

Adding the two equations for Cq and Ci in equation ([8]), we get 



Li 



(9) 



KRTf + 1 

Co + Ci equals to the concentration of bound receptors. We can thus assess our assumption 
that the free receptors outnumber greatly the bound receptors. Substituting typical concentrations 
Lf and assuming kR ~ 1, we do get that Cq + Ci <ti R. The condition is stretched a bit for self 
ligands (which are more numerous), but the analytical results derived are nonetheless in excellent 
agreement with numerical results for all but unrealistically large ligand concentrations. 
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Now, we have at steady state 

aKCo 



and K 



uKtCq 



(10) 



h - l+(f)(7o ^ (r7^ + 6)(l+(f)Co). 

Observe how Ci oc KCq and K oc Co~^ at large Cq. This shows right away that Ci will 
be independent of Co and thus of L/ at large enough ligand concentrations. One can substitute 
equation ( fTO] ) in equation (|9]) to obtain 



Co 1 + 



T 



(11) 



(r/i + 6) (1 + (f ) Co) y \i^RTf + l, 
This is a quadratic equation for Co. It has only one positive solution, whose exact form is not 
important, but which reduces for large Lfio 



Co^ 



uRTf 



^uRTf + 1 



C^aKrTf 
- 1 , for large. 



1 + bT 



f 



We defined C* = e(5 ^. For large Lf "we then have, considering equation (10): 



(12) 



l + 6r/ 

(7i monotonically approaches this hmit from below. 

The last important quantity to solve for is the scale of concentration of Lf for which Ci is close 
to its asymptotic value. A natural measure is the ligand concentration for which the output reaches 
one half of its asymptotic value. We see from equation ( [T0| ) that this occurs when Cq = C^. 
Solving for Lf in such case yields 

More generally, one can show that L/ a defined through Ci{Lf = I/j a) = A • Ci(Lj ^ oo) 



equals, from equation ( 10) and (11) (and taking 6 = to compare with the parallel sorting case) 



/ kRtj + 1 



(13) 



B. Parallel Adaptive Sorting 

Below is the general case with upstream and downstream kinetic proofreading of the adaptive 
module. is the total number of steps in the signalling cascade, m denotes the complex activating 



the adaptive node. The symbols have the same meaning as in Sec. |II A[ The only restriction is that 
< m < N-2. 
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N 



N 



\ i=0 
<P {Cn-l — Cn 



a 



Dn-1 
Dn 

k 



i=0 

for 1< n < N - 2, 



CN-^ — 



4)Cn-2- {aK + ry^) Cn-1, 

- E + ^^]) - E - + ^«"') 

<P {Dn-i - Dn) - T-^Dn for 1 < n < iV - 2, 
(/)Dn-2- (aK + T-^) Dn-1, 
aKDN-i - t'^Dn, 
-S{C^ + DjK + e{KT-K). 



Do = K 

Dr,=<J) 



(14) 




The schematic representation of this network is shown in Fig. [5} 

free free ^b ^b ^b ^b ^ 

x:^ b b b b b b 



FIG. 5. Schematic representation of parallel adaptive sorting. Foreign complexes (C^, z > 0) decay to R^^^^ 
and with rate r^^ (not shown). Self complexes (Di, i > 0) decay to i?^^^^ and L^^^^ with rate (not 
shown). The output here is Cn Dn (circled). All arrows with no specified enzyme are unregulated (fixed 
rate). 

The only new variable is 0, which is the (unregulated) phosphorylation rate down the cascade. 
In order to make analytic progress, we take the dephosphorylation rate 6 = 0. The qualitative 



features of the solutions do not depend on that assumption (see Sec. |IIC| for a discussion). 

Observe that the only form of coupling between the self and foreign is through receptor se- 
questration (the term R — + ^^]) phosphorylation of the kinase (shown as dashed 



arrows in Fig. l5 ). As before, we take R^^^^ = R — J2iLo [Ci + A] ^ in the above equations. 
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We consider uniquely the steady state. Summing the equations for the Ci, one obtains 



1=0 

Now, note that 



= ( .J^IL_ ) c^_^ for 1 < n < - 2. 



To alleviate the notation slightly, define 7/ = (j)Tf{(j)Tf + One can then repeat the above 
to obtain C„ in terms of Co, that is 

C„ = 7^Co for 1 < n < iV - 2. (16) 



One can also easily obtain that from the original system of equations (14) 



and f K \ 

Cn = aKr^C, = (^J^j </>r/7r^Co. (18) 

We can substitute equations ( [T6| ), ( [T7| ) and ( fTS] ) into equation ([15]). Interestingly, the factor of 
K (containing the coupling between the two types of ligands) cancels from this equation. We have 

From the definition of 7, we can further simplify this to 

C. . (^\ . ,19) 



The above derivation did not depend the fact that we were specifically considering foreign 
ligands. The above expression is then also valid for (with Lf ^ Lg and r/ r^). 

We are interested in the output Cn + D^. To solve for it, we must determine the steady-state 
value of K. From our system of equation, we read 

K= r = VI — ^- (20) 

l + c:;\Cm + D,r^) 1 + c*"' (t^ c^o + T^^o) 

As before, = e5~^. Since and Co are both known functions of and Lf respectively. 



then so is K. We can rewrite equation ( 18 ) as 



This is equation (7) from the paper. The equation for the same (foreign parameters 

changed to self ones). and Co are known, so the above represent a closed solution. 

Suppose we consider only one type of ligands (e.g. foreign). Then, since Cq (x Lf, we have 
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for large Lf that Cn tends to 

/I \ N-2-m 

Cn ^ {aKrC^c^) r] [j^^^j ^/ ^^^g^- (22) 

Observe that Cat is monotonic in Cq. The above thus constitute an upper bound on the output. 

We see then that Cn is at least quadratic in r . It follows that in the presence of self ligands (L^ > 0) 

the contribution of output from self ligands is then truly negligible, since <C r/. 

It is desirable to determine the scale of at which this asymptotic limit is attained. We can 

compute the concentration in foreign ligands, denoted Lf{A), required to attain a fraction A of the 

asymptotic output. One can arrive from equation ( [T9| ) and ( [21] ) at 

Notice that the result does not depend on N, the total number of steps in the cascade. Since 
we want a small Lf{A), (j) cannot be arbitrarily small or large as Lf{A) grows as ^ oc and as 

The contribution to Lf{A) due to > will be small given that 7s is small (see Sec 



IV 



for a discussion). 7s is a monotonically increasing function of 0, with value TsT^^ for = 
and tending to 1 for ^ oc. At = r^^, we have that Jslj'^ < "^^s^f^ ^ 1 already. 

The intrinsic contribution (that present if Lg = 0) has a non-monotonic behaviour in 0. If we 
plausibly assume aKx = (taking all steps in cascade to have same default rate of phospho- 
rylation), then, for a given C*, the minimum of Lf{A) occurs for cj) = \mTj^. Taking 
could decrease self antagonism by a factor of two. However, this would be at the detriment of the 
intrinsic sensitivity, which would end up dominating at very small (j). cj) ^ ry^ sets the appropriate 
scale for a maximal sensitivity (optimizing intrinsic sensitivity and mitigating self antagonism). 



C. Effects of Non-Zero Cascade Dephosphorylation Rate 

The analytic calculation for the parallel adaptive sorting module was performed assuming a 
vanishing dephosphorylation rate (denoted by b) down the signalling cascade. In order to better 
understand the effect of 6 > 0, we can first consider the simple adaptive sorting module. This case 
can be solved with b > 0. From equation ( [T2| ), we have 



0(r/,6) _ frA fl + br, 



whereO(r, 6) = lim Ci(L,r, 6). 



We see that the ratio in output arising from foreign and critical ligands is a monotonically 
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decreasing function of b. This is intuitively clear: in the limit of large 6, the dominant opposing 
force to the production of the output is no longer unbinding of the ligand-receptor complex (which 
depends on the parameter of interest, r), but rather the backward reaction rate down the cascade. 
It is thus expected that the output concentration strongly depends on b and only weakly on r in the 
limit where brf ^ 1. 

These observations for simple adaptive sorting actually apply to the parallel sorting module. 
For 6r/ :» 1, the output due to different binding time are indistinguishable and all specificity is 
lost. Figure |6] below illustrates the results of numerical integration for a chosen set of parameters. 
Qualitatively, the behaviour is the same as simple adaptive sorting. 
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10"^ 10^' 
Dephosphorylation Rate (b) 



FIG. 6. Left panel shows output versus ligand for Tf and Tc (full and dashed lines respectively) for 
different values of the dephosphorylation rate b. Right panel shows the ratio in foreign and critical output 
(large ligand concentration asymptotic value) as a function of dephosphorylation rate evaluated numerically 
(solid line), b values from left panel are identified as dots. The dashed line is the 6 = limit (from equation 



(|22|)). Note how dephosphorylation degrades specificity: at large dephosphorylation rate down the cascade, 
the output from foreign and critical ligands tend to the same value. Parameters: = 4, m = 2, = 0, 



i? = 3 X 10^, K = 10- 



0.3, Kr = 1,5 = land 6 = 0.5. 



As discussed in Sec. [V} assuming a low dephosphorylation rate for bound pMHC-TCR and 
a high dephosphorylation rate for unbound pMHC-TCR is consistent with the actual biological 
system of interest even there might appear to exist tension between these requirements. 



III. DISCUSSION OF THE SPECIFICITY/SENSITIVITY TRADE-OFF (INTRINSIC) 



As we detail below, the biochemical networks considered in the previous section cannot be 
arbitrarily sensitive (respond at very low concentration of ligand) and be specific (have very differ- 
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ent output concentrations for not very different ligands). The limits derived concern the "intrinsic" 
properties of the networks (i.e. in absence of self ligands). A discussion of the effects of self 
ligands on sensitivity is presented in the next section. 



A. Simple Adaptive Sorting 

The threshold in output at which the system responds is an important variable of our model. 
To derive the minimal possible sensitivity, it is desirable to take the threshold to be as small as 
possible while still retaining the self-consistency of our framework. 

Critical agonist ligands are defined to be the ligands for which an arbitrarily high concentration 
of ligands does not trigger response. One must then demand that the threshold be above the 
maximum output value that can be attained by the critical ligands (r^ ~ 3 s): 

Threshold > Asymptotic output critical critical ligands = — — - — - . 

1 + bTc 

This is a hard bound needed for the internal consistency of our model. If the threshold were 
placed lower, it would imply that critical ligands could trigger a false positive. Note that in practice, 
in the presence of fluctuations, the threshold must be placed higher to avoid having fluctuations 



cause response for critical ligands. This point is treated in detail in Sec. VL Here, we are concerned 



with strict lower bounds on the sensitivity and in particular in Sec. [IV] on the minimal effect of self 
ligands on this sensitivity. 



We can use equation ( [T3| ) to compute the foreign ligand concentration needed to reach that 
minimal threshold. To do so, we must take (taking 6 = gives us a strict minimum) 

^ _ Ci{Lc oo) _ Tc 
Ci{Lf^oo) Tf' 
Substituting this A in equation ([13]) yields 

L^^n = a (^^1^) { 7^ + «^Tre} (24) 

Lmin Specifies the minimum sensitivity of the model: ligands with r = Tf will not trigger 
response at lower concentrations. L^in is optimally as small as possible, since we would like the 
response to be triggered by very few foreign ligands. 

Concerning specificity, note that from equation ( [T2] ) we must have 

"^''^i ^ < CaKr. (25) 
dr 

The above is a measure in our model of the specificity, i.e., of how fast the asymptotic output 
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concentration varies with the binding time. Optimally, we would like this quantity as large as 
possible. Small difference in ligand types (characterised by r) could then be mapped to large 
differences in output concentrations. 



Looking at equations ( |24| ) and ( [25] ), we see however that the requirements of arbitrarily low 
Lmin and high dCi/dr cannot be met. In effect, we have the constraint (eliminating from the 
equations) 

The prefactor in front of dCi/dr quantifies the intrinsic trade-off between specificity and sen- 
sitivity. 



B. Parallel Adaptive Sorting 

It is possible to apply the same idea in the case of adaptive sorting with upstream and down- 



stream kinetic proofreading (Sec. II B ). The notion of minimal threshold presented in the previous 



sub-section holds (i.e., that the threshold must be larger than the asymptotic output concentration 



corresponding to critical agonist ligands). In particular, using equation ([23]) with the appropriate 
A, one can obtain Lj^in in the more general case. 

As before, it is possible to eliminate one parameter (C*) in favour of dCN/dr in the expression 
for Ljnin to obtain a bound of the form 

Lmin>^N,m{(t>,C^KT) ( , whcrC 



dr 



nr / o —9 \ / —1 , \ N — m — 2 



Here, tm = 0.5(r/ + Tc). Introducing tm is needed since in the general case dC^/dr depends 
on r . To get a representative value of dCN / dr, we take r to be the mean of r/ and Tc. F Ar,m (defined 
only for m + 2 < A^) quantifies the severity of the trade-off between the intrinsic specificity 
and sensitivity. Clearly, the specific form of the function is not critical. What is important is 
that at a fixed m, increasing A^ decreases Tn^^' Similarly, at a fixed A^, increasing m increases 
r7v,m- Precisely, one can show that Tn^^ = (1 + 0"^^c"^)^ ^n,o > ^n^o- Such behaviour can be 
traced back to the form of the solution in the general case, equation ( [22] ), which contains a factor 
{(/)r{l + (/)r)"^}^"^"^. Hence, larger N — m implies stronger dependency on r, making it easier 
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to evade the intrinsic specificity- sensitivity trade-off. Fig. [7] illustrates these trends. In the figure, 
we take aKx = to reduce the number of parameters. This simply means that the phosphorylation 
catalyzed by kinase K is at the same default rate than the other ones in the cascade. To be clear: 
small r implies less trade-off (better). 




in - 2 



FIG. 7. Top panel shows Fg^o as a function of (p. We explicitly assume aKx = 0. Lower panels show 
rN,m normalized by Fs^o- From left to right: m = 0, 1 and 2. Observe how, at fixed m, FAr,m decreases 
with increasing N. Moreover, notice that at fixed N, FAr,m increases with m. nR = 3 in the above. 



The m = case is important as it shows how > 1 drastically relaxes the trade-off. For 
m > 0, the improvement is still substantial. 



IV. EFFECT OF SELF LIGANDS ON SENSITIVITY 



As heuristically shown in the paper, self ligands have deleterious effects on the sensitivity of 
the adaptive sorting networks. Fig. [8] compares heuristic conditions (5) and (7) from the paper 
to numerical solutions for L1/2 (foreign ligand concentration required to reach half maximum of 
asymptotic output). 

In the following sub-sections, we derive mathematically rigorous bounds for the loss in sensitiv- 
ity due to self ligands in both the simple and parallel adaptive sorting modules. This is essentially 
an elaboration of the results of equation (5) and (7) from the paper. 
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FIG. 8. Ligand concentration required to reach half of maximum output concentration. Results for both 
simple adaptive sorting (green) and parallel adaptive sorting (purple, example with N = 3, m = 1) are 
shown. Parameters (same for both simple and parallel adaptive sorting): i? = 3 x 10^, n = 10~^, = « = 
0.25, 6 = 0, Kt = 1, 5 = 1 and e = 0.5. Dots are obtained from numerical integration of the network 
equations. Lines are the quantities appearing in equations (5) (green) and (7) (purple) in the paper. Note 



that at low Lg values, the "intrinsic" sensitivity dominates (see Sec. III). 



A. Simple Adaptive Sorting 

In the presence of two types of ligands, trying to solve for the steady-state of equation ([8]) yields 
a quartic equation. The exact solution yields no insight. Instead, we derive a strict lower bounds 



on the shift in the minimum ligand value required for response (defined in Sec. [HI]) due to the 
presence of self ligands. 

In order to obtain a lower bound on the loss of sensitivity due to self ligands, we first obtain an 
upper bound on Ci. To do so we need a lower bound on Dq. This is because Dq appears in the 
equation for Ci (through the kinase which couples the two types of ligands). 

Dq is a monotonically decreasing function of the kinase concentration. Setting K to its maxi- 
mum value Kt yields 



\K.RTs + 1 J (-, _L_ aKT ^ \K.RTs + 1 J (1 + aKrTs 

The lower bound in Dq then gives us an upper bound on the value of K in the steady-state 
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1 + C^^ (Co + Dq^Iow) 

This finally gives us an upper bound on the concentration of Ci, satisfied for all concentrations 
of self ligands 

aKCo Kt I o^Cq 



' Tj^ + h- 1 + (Co + D^^io^) \tJ^ + 
We now see that our bound on Ci is a monotonically increasing function of Cq. It is possible 

to bound Co above . ^ 

C,<C, + C,= (;J^ j ^/ = ^o,«p. 

This gives us our final rigorous upper bound on Ci 



C, < ^rr^ , ( 1 = C, 



up- 



1 + C* {Co^up + Dq^iow) +b 
Note that the presence of self ligand will contribute through Di to the output concentration 

(taken to be Ci + Di). Because of the smallness of t^^Ts, this contribution is however negligible. 

We can obtain via Ci^up the minimum foreign ligand concentration needed to reach the minimal 
threshold value in the presence of self ligands. This allows us to obtain a lower bound on the loss 
of sensitivity in presence of self ligands ALj^iniLs) = Lj^iniLs) — L^niniLg = 0). We get 

Tc 1 / 1 + tvRrf \ f Lg 



Ai™...(i.) > [fj (26) 

The term in {} depends on the specific choice of threshold. The second term is slightly smaller 
than the scale stated in the paper (equation (5)) because of aKTTg (which was neglected in our 
heuristic argument). Numerical integration show that this bound is indeed satisfied (see Fig. [9]). It 
is known that ^ 10^ on the surface of antigen presenting cells. Hence, our simple but rigorous 
bounds tells us right away that if the adaptive node is activated by Co, then ALj^in^Lg) > 2000 
where data has response around L ^ 5. This means that in the presence of multiple self ligands, 
the ability of this simple model to respond at low foreign ligand concentrations is annihilated. 

B. Parallel Adaptive Sorting 



Interestingly, the more complicated case (Sec. II B) is solvable exactly. In particular, equa 



tion (23 ) allows us to compute ALmin{Ls) defined in the previous sub-section. Indeed, ALmin{Ls) 



is the contribution coming from the self ligands (Dq) in equation (23), so 
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Substituting for the expression for Dq in terms of Lg (equation ([19])), and using A consistent 
with the minimal possible threshold, one arrives at (neglecting terms of order since A ~ 0.1) 




N-m-2 




m+1 



(27) 



The term in {} is specific to our choice of threshold, but the second term corresponds to the 
scale given in the paper (equation (7)). The relevant term is the last factor before Lg. Since 
(j) ^ tJ^, this is a small parameter. We see thus that increasing m (i.e. moving the activation 
of the adaptive complex downstream) strongly decreases the ability of the non-agonist ligand to 
desensitize response. Comparison of equation ([27]) to numerical evaluation is illustrated in Fig. [9| 
below. 
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FIG. 9. Left panel shows schematically the definition of ALmin as a quantifier of loss of sensitivity due to 
self ligands. Full lines are for = and the dashed line for Lg > 0. Right panel shows numerical calcu- 



lation of AL together with analytical results of equations ( [261 ) (green) and ( |27| ) (purple) for simple adaptive 
sorting and parallel adapative sorting respectively (calculations with = 3 and m = 1). Parameters are 
the same as in Fig. [8j There is slight disagreement at large Lg due to receptor saturation. 

It is not possible here to quantify a trade-off between specificity and sensitivity as cleanly as 
for the intrinsic case (Sec. [III]). Such a trade-off between minimizing the self antagonism and 
maximizing specificity is nevertheless present. To see this, note that the only parameter that we 



could change in equation (27) is (f). ALrnin{Lg) is a monotonically decreasing function of 0. 



Our result for Cn (equation ([22|)) show that the output specificity increases with increasing 0. As 
before, the two requirements of decreasing AL^^^(Ls) and increasing specificity are incompatible, 
but increasing for fixed m relaxes the trade-off. 



26 



V. CONSEQUENCES OF FINITE UNBOUND RECEPTOR DEPHOSPHORYLATION RATE 



An important assumption of our model is that the dephosphorylation of the internal section of 
unbound TCR is fast. This leads to all complexes in the signalling cascade decaying to R^^^^ and 
L^^^^ with a rate equal to the inverse binding time (e.g. Fig. 2 (d) of the paper). How specifically 
our model depends on this assumption needs to be assessed. 

To this end, consider the conservative model where dephosphorylation of the internal section 
of the unbound TCR occurs distributively (worst case scenario) with finite rate /5. Assuming that 



the rest of the biochemistry is unaffected, we have the network illustrated in Fig. [TO 

-1 



Ro+L 



free 



J_ b b b b /-X "I 

T^'Co^^ C.:^-:^ I Bound TCR 

P J ^ ... ^ ^ \ Unbound TCR 



FIG. 10. Parallel adaptive sorting scheme with finite dephosphorylation /3 of unbound TCR. Rn denotes 
the unbound TCR phosphorylated n times. As in the parallel adaptive sorting module (Fig. [5]), kinase K is 
deactivated by Cm and Rm (not shown for clarity). Only one type of ligand is shown. 

It is intuitively clear that for ;5 oc, the new model reduces to the parallel adaptive sorting 



scheme described in Sec. |IIB[ We consider the situation /5 < oc to understand limitations of our 
description. 

The first thing to notice with /3 < oc is that even in the absence of pMHC molecules, there 
will be basal levels of modified TCR (i.e., Rn > for all n). In the absence of ligands, it is 
straightforward to show that 



1 - 4,13-^ 



< 1 (for n > 0). 



(28) 



In the presence of ligands, the steady-state value of Rn is more complicated although we expect 
the above expression to be valid at low ligand concentrations. We are interested in determining 
whether or not the properties of the network is strongly compromised by /3 < oc. One property of 
particular importance is the fact that the network must be sensitive to low concentrations of foreign 
ligands. In that regard, note that the equation for the steady- state value of the kinase K becomes 
in the present case, using equation ([28]) (cf. equation ([20|)) 



27 



K 



T 



Krj 



1 + C*-' {Cm + Rm) 1 + C^' {Cm + 0-/5 — i?) ' 

We see that just as introducing another type of hgand (self hgands, see equation (4) of paper), 
introducing a finite /3 leads to a shift in the ligand concentration required to reach the adaptive 
regime. Crudely, taking our low ligand expression for i?^, we expect that a finite /3 will lead 
to a shift in the half-maximum ligand concentration given by (neglecting a term of order 1 since 
hiRr > 1 for r > 1 s) 



Cm ^ R 



ALi/2(/3) 



/3r 



(29) 



Figure [TT] illustrate this effect for a given set of parameters. Equation ([29]) captures quite well 
the effect of a finite /3 on the sensitivity of the model. The larger the smaller the decrease in 
sensitivity. Interestingly, we see from equation ([29]) that increasing m (i.e., moving the position of 
the activation of the adaptive sorting module downstream) helps to reduce the "basal" deactivation 
of the kinase K by unbound TCR. See below for an opposing effect however. 
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FIG. 11. Left panel shows output versus ligand relationship for different values of /3 (color coded to right 
panel). Inset shows the "basal" output level as a function of 13. Right panel shows the shift in L1/2 due to /3 



(shown in left panel). Points come from numerical integration and dashed line is equation ( [29| ). Parameters: 

= 4, m = 2, T = T/ = 10 s, = 0, = 3 X 10^, = 10"^, = a = 0.3, 6 = 0, = 1, 5 = 1 and 
6 = 0.5. 

An interesting feature of the network with /5 < oc is that at low ligand concentration, the 
output shows a non-monotonic behaviour in /3 (Fig.[TTJ left panel inset). One way to conceptually 
think about this is that with finite /3, the unbound receptor behave essentially as another ligand 



type with binding time ^ f3 ^. From equation (29), we see that for /3 > cj), the main effect of 
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modified unbound receptor is to desensitize response. This explains the initial decrease in output 



as /3 decreases. As /3 gets on the same order of magnitude as 0, we see from equation ( [28] ) that the 
unbound receptors will actually begin to contribute to the output itself. This explains the increase 
in output concentration for low /3. 

For our model to remain sensitive in the presence of finite unbound dephosphorylation, we 



must demand that the unbound TCR do not deactivate the response. From equation (29), this 



corresponds to ALi/2{(3) <C 1, or to ;5 > (0 + r"^) ^^^(rT^r). 

It is also important to mention that the potency of self ligands to desensitize response will 
be enhanced by a finite /3. This is qualitatively easy to understand: a finite /3 leads to more 
intermediate complexes in the cascade as a result of a slower rate to R^^^^ and L^^^^. In particular, 
a finite /3 will increase L^^, the complex arising from self ligands responsible for the deactivation 
of the kinase. This means that the condition Cm ^ required to reach the asymptotic output 
concentration will be harder to achieve, implying a lower sensitivity to foreign ligands. This effect 
arises as Cm is not as strongly affected as Dm by /3 < oc because Tf Tg. To see this, one 
can very roughly approximate the effect of a finite /3 as an increase in the binding time of the 
complexes in the cascade. To return to R^^^^ and L^^^^, complex n must not only unbind, but also 
undergo n dephosphorylations. So, neglecting rebinding possibilities (reasonable for large one 
has effectively that r ^ r + n/3~^ for complex n, where /3~^ is the mean time needed for one 
dephosphorylation of an unbound receptor. It follows that the relative effective change in r will 
be more important the smaller r is, implying a more pronounced effect of /3 < oc for complexes 
arising from self ligands. From this simple estimate, we know that this enhanced desensitization 



due to self ligands should be small provided f3 :» mr^ ^. Numerical simulations (e.g. Fig. 12) 



confirm this. This restriction for the conclusions of our model to be valid is in addition to the 



previous one arising from equation ([29|)). 

Note that from Sec. |IIC[ 6 <C r was also required for optimal specificity. There thus seems 
to be a tension between these two requirements: small dephosphorylation rate of unbound TCR 
but large dephosphorylation rate of bound TCR. However, it must be stressed that this model 
is fully consistent with the kinetic segregation model ifTTIl , which has recently received strong 
experimental support ll2T1l . In that model, the binding of the pMHC to the TCR approaches the 
cells' membranes, thereby pushing phosphatases with large extracellular domain away from the 
region of binding. As the pMHC-TCR complex unbinds, the phosphatase can have access to the 
TCR again. We thus expect, in a very crude sense, that /3 be large and b be small in a model of early 
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FIG. 12. Shown is the loss of sensitivity (quantified through L1/2) to foreign output as a function of 



^ in the presence (red) and absence of self (magenta, essentially the same data as in Fig. [TT]). The blue 
curve shows the behaviour of L1/2 arising from self only (subtraction of L1/2 with L5 = from that with 
Lg = 10^). We see that as f3 decreases, the potency of self ligands to desensitize (larger L1/2) response 
increases. Parameters as in Fig.[TT] except for Lg. 



immune signal transduction incorporating kinetic segregation. Complex membrane biophysics 
would need to be taken into account to do full justice to the kinetic segregation model. This 
is beyond the scope of the present paper. Our model is in a way a zeroth order approximation 
consistent with the essence of kinetic segregation. 



VI. ROBUSTNESS OF ADAPTIVE SORTING TO STOCHASTICITY 

Our discussion thus far was purely in deterministic terms. Given the low concentrations of 
ligands considered, such perspective is incomplete. 

It is crucial to verify that the performance of adaptive sorting is not compromised by stochas- 
ticity. In what follows, we perform all our stochastic simulation of the adaptive sorting networks 
with the Gillespie algorithm ll22l . 



A. General Considerations 



In the present work, we are specifically concerned with the steady- state output concentration 
0(r, L) of the signalling pathway, as a function of the ligand concentration and binding time. 
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In order to make progress in the understanding of the stochastic properties of the displayed 
solutions, we make the following reasonable assumption for a signalling cascade: The output con- 
centration is lower than the concentration of other species. If this is true, then it is possible to treat 
relevant effects of fluctuations by focusing specifically on the output. A natural way to do this is 
to assume that the output follows a Poisson birth-death process, with production (p) and degra- 
dation (S) rates given by the corresponding deterministic rates. Doing this neglects correlation 
between the output and other species, but should be valid in the limit that our initial assumption 
(output concentration smaller than other species' concentration) is valid. This approximation is 



schematically illustrated in Fig. 13 for the specific example of a parallel sorting scheme. 




ftts free ^ ^b ^b 
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FIG. 13. We approximate the full network as a simple Poisson birth-death process for the output. This 
approximation should be valid provided the concentration of other molecules is large compared to the 
output. The subscript "ss" denotes steady-state values. denotes the fictitious pool of molecules from 
which Cn is produced and degraded to. 



We can further approximate our Poisson birth-death process by a simple Langevin scheme Il23]| . 
This approximation is valid at large molecular number, which is not strictly the regime of interest. 
We need this to make analytical progress. The stochastic differential equation is 



^ = -SO + V{t)\l2p + 5d ^ -SO + r(t)^2^, where = + pS-\ 

T{t) is Gaussian white noise with unit variance. The initial condition is that 0{t = 0) = —pS~^. 
By construction, we have that pS~^ = O^^^(L), the deterministic steady-state output concentration. 
It is important to stress that in our scheme this is independent of the specific topology or param- 
eters of the underlying network. The second approximate equality above comes from neglecting 
O in the square root. Since in the steady-state, (O) = 0, this amounts to doing a first order 
approximation. (•) denotes an ensemble average. 

In the end, we see that the full network can be approximated by a Ornstein-Uhlenbeck process, 
whose stochastic properties are well known. In particular, note that |[24ll 



{0{t)) = -0'^\L)e 



-St 



(30) 
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= {d{t)){d{t')) + 0^'\L)e-^^'-''^ (l - e-^^''^ {t > t'). (31) 

Since we are considering low number of output molecules, the output's time course will be 
highly fluctuating on short time scales. As a result, it is not possible to make decisions concerning 
the nature of the bound ligands (critical or foreign) given observations of the output alone. 

Motivated by the biological system, we introduce species A downstream of the output respon- 
sible for time averaging the output ifTTll . In our model, the thresholding mechanism acts on A, not 
on O. It is reasonable to assume that the discrimination and amplification steps are realised by 
different modules in the actual biological system. We take 

i = AO - T-^A, 

where T is effectively the averaging time. We suppose in what follows that the kinetics of A is 

slow compared to that of the output {5T ^ 1). A is a kinetic first order constant that sets the 

scale of concentration A. Many molecular species could play the role of A in signal transduction 

of the immune system and most of them have typical concentrations above ^10^ [7J. It is then 

reasonable to assume that A is large enough that stochasticity of A will be entirely dominated by 

the stochasticity of the output O. We therefore integrate for A deterministically. With ^4(0, L) = 0, 

we can formally write the solution of A{t, L) (note that A depends on the ligand concentration) as 

function of 0(t, L), . 

A{t,L)=K e-^'-''^'^0{t',L)dt'. (32) 

^0 



With our approximations, equations (30), (31 ) and (32) together with our definition of O allow 



us to compute the stochastic properties of A{t,L). One finds (recaUing that p5 ^ = 0^^^(r, L), the 
deterministic steady- state output concentration) 

{A{t, L)) = 0^'\L)AT{1 - e-'l^ + fe-'l^ (l - e''!^) } where f^T[bT- 1)"' , (33) 

= L)f) - {A{t^ L)f ^ 0'-{L) (1 - e-^^/^) . (34) 

We neglected transient terms smaller by a factor of 5~^T~^ in our expression for a\^^ Given 
the level of approximation we are after, the exact form of the variance of A{t, L) is not more 
useful than the above compact expression. From these results, we can see for what value of T 
discrimination between critical and foreign ligands become possible with parallel adaptive sorting. 

To do so, we must first determine a plausible threshold value. Optimally, the threshold should 
be positioned at a value as low as possible (to be sensitive to as little foreign ligands as possible), 
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but not sufficiently low as to trigger response with critical ligands (no false positive responses). 
Since the output increases monotonically with ligand concentration in an adaptive sorting scheme, 
we can consider the large Lc limit for 0^^\Lc), denoted by 0^^\ We can set our threshold 9 to be 



eiij) = (A(r,)) + ijaA^ = AT | Of + ^^/^^ | , (35) 

where quantifies how conservative the system needs to be about auto-immunity. Ac above de- 
notes the large time and ligand concentration limit of A for critical ligands. Larger ^Jj implies less 
frequent false positives. The underlying biochemical network will provide a suitable discrimina- 
tion mechanism even in the presence of noise provided one can still find 

(^A{Lf) 

on the order of one or larger for even very small foreign ligand concentrations. Z{Lf) tells us 
how far from the threshold (measured in units of typical fluctuations of A) the output arising from 



foreign ligands is at concentration of Lj. Indeed, from equations (33), (34) and the definition of 



the threshold 9, equation (35), one gets 



Z(i, L,) « ^0^(L,) (1 + Ti,) {l - ^} - iy (36) 

It is approximate only as it neglects terms of order bT~^. We can evaluate this quantity knowing 
the model for the signalling cascade. It is worthwhile to substitute numbers. In our model, ^ is 
essentially the unbinding rate r"^. Taking O^^^ = 0.5 and Of\Lf = 5) = 1 are conservative 



estimates within our models. With r/ = 10 s and Tc = 3 s, we have Z ^ 0.5^1 + O.IT — 0.4?/;. 
For T = 60 s and = 3, we already have Z ^ 0.1, meaning that we expect more than half of the 
cells to respond in a time- scale of one minute when exposed to on the order of 5 foreign ligands. 



The important point to realize in equation (36) is that as T increases, the relative importance 



of (encoding the conservatism of the system against auto-immune reaction) decreases. This 
simply makes explicit the fact that relative fluctuations decrease as T"^/^. One cannot take T to be 
arbitrarily large, as this would lead to response time too long to be useful in a biological system. 
However, numerical simulations (next sections) show that T ^ 60 s is an appropriate compromise 
between speed and time averaging of fluctuations. Interestingly, this is precisely the time scale at 
which cells of the immune system respond in the presence of foreign ligands f7l|. 

One might object that for low concentrations, the output will fluctuate above and below the 
threshold. This could in principle incapacitate the downstream thresholding mechanism. These 
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fluctuations occur on time scales that are slow (T, the relaxation time of ^4). There is good em- 
pirical evidence that once the molecular equivalent of A is exceeds a threshold in T-cells, further 
molecular reactions actually deactivate the incoherent feed-forward loop [7 J (i.e., reactivating the 
kinase K in our model). This in turn leads to a sharp increase in the output, making fluctuations 
irrelevant. It is reasonable that the kinetics of this additional reaction are faster than that of A. As a 
result, this means that response would be triggered provided the output stayed above threshold for 
a sufficiently long time, dictated by the kinetics of the mechanism downstream of A. An approxi- 
mation is then to take response to occur as soon as the output exceeds threshold. This is the strategy 
used to generate Fig. 4 (d) of the paper. Observe that such strategy makes the system potentially 
vulnerable to critical ligands. In spite of this, very few cells 1% in our time window) respond 
to critical ligands, even at large concentrations (see Fig. 4 (d) of the paper). In the examples below, 
we instead display the fraction of cells above threshold at a given time (i.e. not assuming anything 
about the time-scale of the thresholding mechanism). The qualitative behaviour of the fraction of 
cells above threshold is very similar to the fraction of cells responding just as threshold is crossed. 
Note that the downstream mechanism is not explicitly taken into account in our model for reasons 
of simplicity. 

Many approximations were needed to arrive at our conclusions. In the following sections, 
we verify numerically for the parallel adaptive module that the conclusions hold. Note that for 
the module activated by the phosphatase (network from Fig. 4(b) in the paper), an analysis of 
noise was performed in ifTTl . It was seen that an adaptive sorting like module (but with a different 
network topology) could have good discriminatory abilities even at low molecular numbers. These 
results indicate that the conclusions drawn from our deterministic results do indeed remain valid 
even in the presence of fluctuations. 



B. Specific Example - Parallel Adaptive Sorting 

We here consider the idealized network performing parallel adaptive sorting (network from 
Fig. 4 (a) from the paper), first in the absence of self ligands (Lg = 0). The first approximation of 
importance is that the output concentration is governed by a Poisson process with production and 



degradation rates equal to the deterministic ones. The top of Figure [14] compares the deterministic 
and stochastic mean concentration in the steady-state (we take a unit volume). 

As expected, the stochastic mean output number differs significantly from the deterministic 
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output at very low ligand concentrations. We see however that there is good agreement even for 



L ^ 3. The bottom of Figure [T4| compares the output molecule number distributions to Poisson 
distribution with rates from the deterministic network (note that the lines are not fit). 



2-0.5 




Oulpul 



OulpLLl 



FIG. 14. Top graph compares deterministic (full lines) and stochastic steady-state output concentrations (A 
for Tf, T for Tc). Lower rows display the output number distribution distributions (squares) for the ligand 
concentrations of the top graph. Poisson distributions (full lines) with the analytical deterministic result for 
the mean of Cat (no fit). Observe how the Poisson distribution with the deterministic mean becomes a valid 
approximation as L increases. Even at L = 5, the Poisson distribution captures well the output distribution. 
Parameters: TV = 4, m = 2, = 0, i? = 3 x 10^, k = 10"^, (j) = 0.3, a = 0.0003, 6 = 0, = 1000, 
5 = 1 and e 0.5. 

For the same parameters, we can also verify that the statistical properties of A as derived com- 
pare well to the full stochastic result. Sample trajectories for A as well as {A{t, L)) ± (JA{t,L) as 
given by equations ( [33] ) and ([34]) with deterministic output given by equation ( [2T] ) are shown in 



Fig.Oand 16 



Note that at moderately high concentrations (L = 160, second row in Fig. [16]), the analytic 
expression is very close to the actual numerical result. This is expected as the correlations be- 
tween the output and other species should decrease in importance as the concentration of ligands 
increases. Agreement of steady- state is also quite good at other concentrations in spite of our 
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numerous approximations (there are discrepancies at intermediate ligand concentrations). 

Systematic discrepancies in the transient behaviour is seen for all concentrations: the analytic 
solution increases too fast at early times. This can be understood simply. In our approximation of 
the full network, (shown in Fig. [13]), we neglected the transient of all the species but the output. 
Specifically, we take the rates governing the Poisson process for the output to be constant. Clearly, 
there is a relaxation time (on the order of r) for these rates to actually attain their steady-state 
value. Therefore, the output steady-state will be attained slower in the full network than in our 
approximation. One way to quantify the discrepancy is to plot the time required for half the 
trajectories to be above the threshold. This is shown in Fig.fTTlfor the concentrations and parameter 



values of Fig. 15 and 16 



Notice how the time for response decreases until a plateau. It is expected in our model that 
extremely high concentrations of foreign ligands do not trigger the response faster than moderate 
concentrations (by construction of adaptive sorting). This is actually seen in experiments (Q], 
see Fig. 3 (c)). In our model, what should decrease the response time is a more "potent" foreign 
ligand, with higher binding time. 

We ran the same simulations in the presence of self ligands: Lg = 10^ with = 0.11 s. (In our 
model, this is as "potent" as Lg = 10^ with = 0.05 s. We need this to speed up the Gillespie 
simulations). We see that just as in the Lg = case, the Poisson approximation is excellent even at 
very low concentration. Our analytic result also capture quite well the exact stochastic behaviour. 



We show below Z{Lf) (equation (36)) in the presence and absence of self. Note again that while 
not perfect, our analytical results capture very well the full stochastic behaviour. 
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FIG. 15. Each row represent a different ligand concentration (same concentrations as in Fig. [14]). The left 
panels present representative time evolutions of A{t). Lines of warm color are with r/ and lines of cool 



colors with Tc. Horizontal black line is the threshold chosen from equation (35) with A = 4. It is the same 
for all rows. Averaging time T is taken to be 60 s. Middle panels show the mean and standard deviations 



of the trajectories from numerical results (blue red rj) and equations (33) and (34) (black). Ensemble 
average over 500 trajectories. Right panels show the fraction of trajectories above threshold at a given time 
(showing those for r/, not more than few trajectories with Tc exceeded sporadically threshold). Note that 
the above is a different quantity than that presented in Fig. 4 (d) of the paper. See discussion. Full lines 
are numerical and dashed lines are from equations ( [33] ) and ([34]). Too few trajectories exceeded threshold 
to have proper statistics for L = 1 (note the scale). Observe the relatively short timescale over which the 



fraction increases. Kinetic parameters are as in Fig. [14 
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FIG. 16. Continuation of Fig. [T5] with remaining concentrations. Observe that L = 160, 500 and 1600 
are practically indistinguishable, as we expect from the main idea of adaptive sorting that output should not 
depend on concentration at large ligand concentrations. 
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FIG. 17. Time required for half of trajectories to exceed threshold (for foreign ligands) as a function 
of ligand concentration. Note the systematic discrepancy of our analytic approximation at large ligand 
concentrations. The analytic result however recovers the qualitative behaviour of the full system. Same 
parameters as in Fig. [15] and [16] 




FIG. 18. Z{tp, Lf) from equation (36). Black line for = and red line for Ls = 10 (r^ = 0.11 s). 



Dots are from full stochastic integrations. Parameters from Fig. [14} [T5] and [16] = 4, T = 60 s). Note 
that already at = 4, there should be ^ 0.5 of response in the absence of self. In the presence of self, this 
goes to L J = 9, which is still a very small number. 

From our analytical expression for Z, we can quickly assess the robustness to parameter choice 



for the parallel adaptive sorting. Figure [19] shows contours of Z for different parameters and dif- 
ferent values of foreign ligands. We see that there is a wide range of parameters where substantial 
response occurs (here with T = 60 s). 
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FIG. 19. Lf) from equation (36) for parallel adaptive sorting module with = 4, m = 2, hcR = 3, 

Lg = 10^, Tg = 0.05 s, aKT = (p, Tf = 10 s, Tc = 3 s. Varied parameters are (p (horizontal axis) and 
(vertical axis). Averaging time 60 s. ^ = 4. Different graphs correspond to different foreign ligand 
concentrations. The green contour corresponds to 50% response. Dark blue curve to ^ 16% response and 
dark red to ^ 86% response. 



VII. ANALYSIS OF TRANSIENT BEHAVIOUR 



We mostly focused w^ith the behaviour of the steady-states in our analysis. This requires jus- 
tification. It important to verify w^hether the transient has properties that could potentially com- 
promise our conclusions. We consider the simple adaptive triage module. We are nov^ interested 
in the full time evolution of (v^e assume only one type of ligand) equation ([8]). Exact analysis is 
not tractable. We obtain approximate expressions. First how^ever, S = Co + Ci can be integrated 
exactly as 



(l - e ^/^^) , v^here = ^ 



T 



Is. 



1 + k.Rt ^ ^ 1 + k.Rt 

We now^ consider the equation for the kinase K. Note that in the present model, Ci^ss <S Co,ss 

(ss, steady-state) at moderately large concentration of ligand. It is then a fair approximation (to be 

verified at the end) to take Co ^ S in the equation for K: 



k{t) ^ - + e) K{t) + cKt 
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It is possible but not very useful to formally write down a solution for K{t). What is important to 
realise is that since S increases monotonically, the rate of decay of K{t) increases monotonically. 
For the purpose of determining the behaviour of Ci (the output), we only wish to determine when 
K{t) is close to Kt and when it is close to its (much lower) steady-state value. It is not hard to 
show then that the following qualitatively captures the behaviour of the kinase: 

K{t) ^ K^e-^^^+^--^« + , 

e + oCq^ss 

where 7y = lim^^ooS(t) = /^i?rL(l+/^i?r)"^. This approximation is valid provided (^/yr^ :» 1. 
We can turn to the behaviour of the output, which is our primary concern. Assuming as before 
that Co ^ S simplifies the equation to 

Ci{t) ^ aK{t)J:{t) - Tc^Ciit), where tc, = r(l + bT)-\ 

The first term on the right hand side represents a time varying rate of production for Ci. The 
second term represents a degradation term with a constant half-life. The important point is that 
since K{t) decreases and T>{t) increases, the production rate will be peaked at some intermediate 
time value. Integrating the equation for Ci leads to 

CM = (^^) t + »Kr f *' (e-."«--<'''E(0 - - ± f C^t). 

V e + oCq^ss J Jo \ ^ + ^^o,s5 / ^Ci Jo 

It must be stressed that the integrand of the second term decays for times larger than r^. If we 
are interested in the behaviour of Ci at large times, then this second term is just a constant which 
can be evaluated by the method of steepest descent. The remaining equation is then straightforward 
to solve (using that r] ^ Co,ss, and that Cq^ss ^ f )• One arrives at an approximate result for Ci{t), 

C,{t) ^ ^ (1 - er^) e-*/^-. + ^^^^ (l - e"*/^-!) for t » r^. (37) 



The first term captures the decay of the initial peak concentration due to the time-varying pro- 
duction rate. The second term represents the relaxation to the equilibrium steady-state. It is 
important to restate the assumptions used to derive the above expression: (1) Stity, ^ 1 to obtain 
our approximate expression for K{t) and (2) Ci <C Co at all times to use T^{t) in place of Co(t) in 
our computations. The latter assumption is self-consistent if the peak value of Ci is much smaller 
than the asymptotic value of C^: aKx <C r]S. Within these restrictions, numerical integration 



agrees with the derived result. Examples are shown in Fig. 20 

The final expression for Ci allows us to check the self-consistency of our approximations. We 
see that the larger 5, then the smaller peak value of Ci. This means that our two assumptions are 
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FIG. 20. Shown in red are numerical integration results for Ci{t) for different values of 6 (with fixed 
ratio I to have the same steady state). The black dashed lines are equations (jsTj). Here, ^ 0.003 
and ^ 0.004. Even out of the range of validity (e.g. 6 = 0.001), the approximation still captures 
qualitatively behaviour of the transient. 



compatible. It is clear however that for very small ligand concentrations, 77 will be small and our 
approximations break down. 

One must keep in mind though that the magnitude of Ci is bounded above by L. Hence, the 
maximum value of Ci cannot be significant at low L. This means that even though our approxi- 
mations do not hold in this regime, only the steady-state behaviour is relevant. 

The whole point of this analysis is to verify that the transient of the output does not have a 
significant behaviour. From our calculations, we see that it can indeed be significant depending 
on the choice of parameters. In particular, for small 5 (scale set by aKxT]'^ and r]~^T^^), Ci 
exhibits a peak that can have a large magnitude. This implies that for small 5, not only the steady- 
state is of importance. In particular, if the immune response is triggered as soon as the output 
reaches a prescribed value (instead of time averaging), then the transient is more important than 
the steady-state. 

More importantly, observe that the peak value of Ci does not strongly depend on r . This means 
that the network looses all its fine r discriminatory properties if 5 is small. It follows that for our 
model to work, the kinetics of the kinase must be faster than the kinetics of the phosphorylation 
of the complex cascade. This fast kinetics ensures a fast decay of which in turn implies that no 
significant amount of Ci can accumulate at short times. In such situation, the important behaviour 
is in the steady state. 
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VIII. EXAMPLES OF EVOLVED NETWORKS 



This section presents networks obtained from our simulations. The raw networks are often 
overly complicated. For clarity, the essential reactions are extracted and shown in simplified net- 
works (reactions that can be removed with no effect on the output behaviour are not shown). 



A. Independent Discrimination - Another Example 



We now display another example of the results obtained via simulation, with only one type of 



ligand present. The simplified network is shown in Fig. [2TJ The output versus ligand relationship 
is shown in Fig. [22} 



FIG. 21. Other example of network found. C^'s are understood to decay to R and L with rate r ^. 




10^ 
[Ligand] 



FIG. 22. Output versus ligand relationship (steady-state) for network of Fig. [21] Note the non-monotonicity 
of the output with ligand concentration. 



As for adaptive sorting, Co is the enzyme phosphorylating some kinase {Ki), In contrast to 
adaptive sorting, Ki does not directly phosphorylate Cq. Instead, Ki is required for the formation 
of which is the kinase responsible for the formation of Ci (the output). Hence, more Co 
implies less Ki which in turn implies less which finally leads to less Ci. This is again the 
same negative feed-forward loop, except with a less direct, "buffered" step. Interestingly, as a 
consequence, the adaptation is not perfect and the response (steady- state output concentration) 
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becomes non-monotonic in ligand concentration as shown in Fig.[22j The adaptive sorting module 
is nevertheless clearly present. 



B. Discrimination in Presence of Self Ligands 

Here we show results of evolutionary simulations when discrimination must be achieved in the 
presence of a large quantity of self ligands (spurious, non-specific ligands). 



1. Example 1 - Phosphatase Activation 



An example of working obtained network is schematized in Fig. [23] (same as Fig. 3 (e) in the 
paper). The full network and parameters can be found in Appendix A. 



Pi 

I 

FIG. 23. Obtained network performing parallel sorting. Complexes Ci are understood to decay to R and L 
with rate t~^. 



Observe first that it is not Co which activates the crucial "adaptive reactions" but rather Ci 
and C3 (dashed circles). The activation of the adaptive module is rewired downstream. Both 
these reaction serve to activate a phosphatase instead of de-activating a kinase (it is still a negative 
feed-forward). In fact, this evolved solution is very close to the realistic model for early immune 
response presented in [17J. Note that the output is not a member of the cascade of complexes, 
but is rather activated by the last element of the cascade (here Cq). This is basically equivalent to 
having the output as Cq directly. The steady-state output versus ligand relationship is presented as 
Fig. 3 (d) in the paper. 

Interestingly, keeping all kinetic parameter constant, it is possible to investigate the effect of 
moving all catalytic activity to the first complex in the cascade, C^. The output ligand relationship 
is shown in Fig. [24} where the dashed line is the output concentration in presence of self ligands. 
The output/ligand relationship does change, but without the presence of self (full lines), this net- 
work would still be regarded as achieving proper discrimination. This is not so in presence of self 
ligands. 
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FIG. 24. Output versus ligand relationship (steady- state) for network of Fig. |23| but with catalytic activity 
of complexes moved to Cq. The dashed line is the output concentration for foreign ligands in presence of 
self ligands (L^ = 10^, = 0.05 s). The effect of self ligands should be compared with Fig. 3 (d) from the 
paper. 



As expected from our analytical analysis of the general case (Sec. |IIB[ ), moving all the catalytic 
activity from complexes to Cq leads to a catastrophic decrease in the sensitivity. The self ligands 
completely inhibit response at low foreign ligand concentrations. 



2. Example 2 - Non-specific Kinase Leading to Non-Monotonic Response 

We display another solution has interesting features in spite its imperfection. The simplified 
network is shown in Fig. [25] 

As before, the complex with catalytic activity is Ci (dashed circles). It is interesting to note that 
here, the adaptive kinase {Ki) is non-specific, phosphorylating two complexes in the cascade {Ci 
and C2). Since Ki ^ L~^, we have that [Output] < {KiYL ^ L~^. The output must ultimately 
decrease at large ligand concentration. This is indeed what is seen in the output/ligand relationship 



displayed in Fig.[26| Interestingly, such non-monotonic behaviour in ligand concentration (loss of 
response at high ligand concentration) is seen experimentally in the immune system (71 [T7|| We 
believe it to be a signature of adaptive sorting for systems with enzymes lacking the biochemical 
specificity needed to act on a single step in the cascade. 

Here again, the adaptive sorting module is identifiable, but the non- specificity of the kinase 
leads to a non-monotonic response instead of a perfectly adaptive response. This non-monotonicity 
is sufficient to separate the output concentrations from the two ligand types over a wide range of 
ligand concentration. 
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FIG. 25. Obtained network displaying parallel triage properties. Complexes Ci are understood to decay to 
R and L with rate t~^. 



As before, for purposes of illustration, it is possible to keep the kinetic rate constants of that 
network fixed but move the catalytic (adaptive) activity from Ci to Co and see the detrimental 



effect of self-ligands on the sensitivity of response. This is shown in Fig. 27 The presence of a 
large quantity of spurious (non-specific) molecules leads a strong degradation of response if the 
adaptive module is activated too early in the cascade. 




10^ 
[Ligand] 



FIG. 26. Output versus ligand relationship (steady-state) for network of Fig. [21] The dashed line is the 
output for foreign ligands in the presence of self ligands. Notice again the little effect of self ligands in spite 
of their being very numerous (10^). 



46 



10" 



10 



3 

o 



10 



10" 



10" 



10' 




10' 

ILigandl 



10^ 



10^ 



FIG. 27. Output versus ligand relationship (steady-state) for network of Fig. |2T] with catalytic activity 
moved to Cq. The dashed line is the output for foreign ligands in the presence of self ligands. The response 
for foreign ligands in presence of self ligands (dashed red) is now difficultly distinguishable from that of 
critical non-agonist (blue) in absence of self ligands. 



IX. APPENDIX - DETAILS FOR NETWORK OF SEC. IVIIIBll 



We provide the details concerning the network of Sec. VIII B 1 , whose output concentration 



versus ligand concentration appears in the paper's Fig. 3 (d). Fig. 28 shows the actual network as 



obtained from evolutionary simulations. One recovers Fig. [23] if the dashed lines are not shown. 
These reactions do not have to be regulated (i.e. catalyzed by enzyme which are modified by 
presence or absence of ligands) for the network to show its distinctive features. 

Initial concentrations and reactions rates follow in tables below. Note that the on-rate is equal 
to tv = 7.02 X 10~^ in the present network. All complexes Q decay to R^^^^ and L^p^ with rate 
ry^. All complexes Di decay to R^^^^ and L^p^ with rate t~^. 
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FIG. 28. Full network as obtained from evolutionary simulations. The reactions denoted by dashed lines 
were not to be important for the behaviour of the system (i.e. enzymes replaced by unregulated ones without 
effect). All complexes Ci decay to the receptor and ligands with rate r^^ (not shown). Complexes from 
self are not shown for clarity. 



Initial Concentrations 



Chemical Species 


Initial Concentration 


Chemical Species 


Initial Concentration 


Ki 


576 


7" free 


Variable 


K2 


475 


Co 










Ci 





K3 


941 


C2 





Ka 


861 


C3 





K5 


348 


C4 





Pi 


439 


C5 





p* 
^1 










P2 


717 




10^ 


^2 





Do 





P3 


957 


Di 





^3 





D2 





P4 


906 


D3 





^free 


3 X 10^ 


D4 
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Phosphorylations 



Kinase 


Species 


Species Phosphorylated 


Rate (xlO ) 


Ki 


Ci 


C2 


1.32 


Ki 


Di 


D2 


1.32 


K2 


C2 


C3 


0.25 


K2 


D2 


D3 


0.25 


Ci 


Pi 


P* 


8.38 


Di 


Pi 


P* 


8.38 


K3 


C3 


C4 


5.88 


K3 


D3 


Da 


5.88 


K3 


C4 


C5 


4.18 


K3 


Da 


D^ 


4.18 


K3 


C5 


Ce 


7.00 


K3 




Dq 


7.00 


C3 


P2 


P* 


4.55 


D3 


P2 


P* 


4.55 


K2 


Co 


Ci 


9.11 


K2 


Do 


Di 


9.11 


Cg 


P3 


^3 


9.40 


Dq 


P3 


p* 


9.40 


Ka 


K2 




3.57 



Dephosphorylations 



Phosphatase 


species 


Species Dephosphorylated 


RateCxlO-"^) 


Pi 


Ki 


K2 


2.77 


^3 


Ci 


Co 


4.24 


P3 


Di 


Do 


4.24 


p. 


p* 
^3 


P3 


8.11 


P3 


^1 


Pi 


1.23 


^1 


C2 


Ci 


8.43 


^1 


D2 


Di 


8.43 


JD* 
^2 


Ca 


C3 


8.73 


^2 


Da 


D3 


8.73 


^2 


Cs 


Ca 


3.34 


p* 
^2 


D^ 


Da 


3.34 


Pi 


Cq 


C5 


2.61 


p* 


Dq 


D5 


2.61 


p* 
^3 


C3 


C2 


4.51 


^3 


D3 


D2 


4.51 


p* 


^2 


P2 


1.21 
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